HIGH-ORDER WAVE PROPAGATION ALGORITHMS FOR HYPERBOLIC SYSTEMS 



DAVID I. KETCHESON*, MATTEO PARSANit, AND RANDALL J. LEVEQUE* 

Abstract. We present a finite volume method that is applicable to hyperbolic PDEs including spatially varying and 
semilinear nonconservative systems. The spatial discretization, like that of the well-known Clawpack software, is based on 
solving Riemann problems and calculating fluctuations (not fluxes). The implementation employs weighted essentially non- 
oscillatory reconstruction in space and strong stability preserving Runge-Kutta integration in time. The method can be extended 
to arbitrarily high order of accuracy and allows a well-balanced implementation for capturing solutions of balance laws near 
steady state. This well-balancing is achieved through the /-wave Riemann solver and a novel wave-slope WENO reconstruction 
procedure. The wide applicability and advantageous properties of the method are demonstrated through numerical examples, 
including problems in nonconservative form, problems with spatially varying fluxes, and problems involving near-equilibrium 
solutions of balance laws. 

1. Introduction. Many important physical phenomena are governed by hyperboHc systems of conser- 
vation laws. In one space dimension the standard conservation law has the form 

qt + .nqh = 0, (1.1) 

where the components of g € M™ are conserved quantities and the components of / : M™ x M™ — M™ are 
the corresponding fluxes. Very many numerical methods have been developed for the solution of (1.1); some 
of the most successful are the high resolution Godunov-type methods based on the use of Riemann solvers 
and nonlinear limiters. These and other methods are generally based on flux- differencing and make explicit 
use of the flux function /. 

Herein we also consider systems with spatially varying flux: 

i^{x)qt + f{q,x)x ^ ip{q,x), (1.2) 

and spatially varying linear systems not in conservation form: 

K{x)qt + A{x)qx ^ ^{q,x), (1.3) 

as well as their two-dimensional extensions. Wave-propagation methods of up to second-order accuracy have 
been developed for such systems in, e.g. [19, 17]. These methods are based on wave-propagation Riemann 
solvers, which compute fluctuations, (i.e. traveling discontinuities) rather than fluxes and thus can be applied 
to (1.2) or (1.3) just as easily as to the conservation law (1.1). 

Second-order methods may often be the best choice in terms of a balance between computational cost 
and desired resolution, especially for problems with solutions dominated by shocks or other discontinuities 
with relatively simple structures between these discontinuities. For problems containing complicated smooth 
solution structures, where the accurate resolution of small scales is required (e.g. simulation of compress- 
ible turbulence, computational aeroacoustics, computational electromagnetism, turbulent combustion etc.), 
schemes with higher order accuracy are desirable. 

The purpose of this work is to present a numerical method that combines the advantages of wave- 
propagation solvers with high order of accuracy. The basic discretization approach was presented already in 
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[13]; here, we give a more detailed presentation and demonstrate the wide applicability of the method. The 
new method combines the notions of wave propagation ([19, 21]) and the method of lines, and can in principle 
be extended to arbitrarily high order accuracy by the use of high order accurate spatial reconstructions and 
high order accurate ordinary differential equation (ODE) solvers. The implementation presented here is 
based on the fifth-order accurate weighted essentially non-oscillatory (WENO) reconstruction and a fourth- 
order accurate strong-stability-preserving Runge-Kutta (RK) scheme. We restrict our attention to problems 
in one or two dimensions, although the method may be extended in a straightforward manner to higher 
dimensions. 

Another approach to high order discretization of hyperbolic PDEs, referred to as the adaptive high 
order derivatives, or ADER method, has been developed by Titarev & Toro [28] and subsequent authors. 
That approach uses the Cauchy-Kovalewski procedure and has the advantage of leading to one-step time 
discretization. The method of lines approach used in the present work seems more straightforward and allows 
manipulation of the method's properties by the use of different time integrators, but requires the evaluation 
of multiple stages per time step. 

A similar class of conservative, well-balanced, and high-order accurate methods has been developed by 
Castro, Gallardo, Pares, and their coauthors; see, e.g. ]2]. Those methods also use WENO reconstruction 
and Runge-Kutta time stepping in conjunction with Riemann solvers, and lead to a discretization with a 
form similar similar to that presented here and in [13]. Those methods have recently been combined with 
the ADER approach; see ]7]. The present method differs in the approach to reconstruction and the kind 
of Riemann solvers used. These differences result in some potential advantages: our method also handles 
systems (1.2)-(1.3) with capacity function k and it can make use of /-wave Riemann solvers ]1] as well as 
wave-slope reconstruction and achieve high order convergence even for some problems with discontinuous 
coefficients. Finally, the method can immediately be applied to very many interesting problems because the 
implementation is based on Clawpack Riemann solvers, which are available for a great variety of hyperbolic 
systems. 

The methods described in this paper are implemented in the software package SharpClaw, which is freely 
available on the web at http://www.clawpack.org. SharpClaw employs the same interface that is used in 
Clawpack [21] for problem specification and setup, as well as for the necessary Riemann solvers. This makes 
it simple to apply SharpClaw to a problem that has been set up in Clawpack. These methods have also been 
incorporated into PyClaw]15], which allows them to be run in parallel on large supercomputers. 

The paper is organized as follows. In Section 2.2, we present Godunov's method for linear hyperbolic 
PDEs in wave propagation form ]21]. This method is extended to high order in Section 2.3 by introducing a 
high-order reconstruction based on cell averages. Generalization to spatially varying nonconservative linear 
systems is presented in 2.4. Further extensions and details of the method are presented in the remainder of 
Section 2. Numerical examples, including application to acoustics, elasticity, and shallow water waves, are 
presented in Section 3. 

2. Semi-discrete wave-propagation. The wave-propagation algorithm was first introduced by LeV- 
eque [19] in 1997 in the framework of high resolution finite volume methods for solving hyperbolic systems 
of equations. The scheme is conservative, second-order accurate in smooth regions, and captures shocks 
without introducing spurious oscillations. In this section, we extend the wave-propagation algorithm to high 
order accuracy through use of high order reconstruction and time marching. For simplicity, we focus on the 
one dimensional scheme and then briefly describe the extension to two dimensions. 

2.1. Riemann Problems and Notation. The notation for Riemann solutions used in this paper 
comes primarily from ]21], and is motivated by consideration of the linear hyperbolic system 



Qt + Aq^ = 0. 
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Fig. 2.1. The wave propagation solution of the Riemann problem. The horizontal axis corresponds to space and the 
vertical axis to time. 



Here q G M™ and A E I]J™>^"\ System (2.1) is said to be hyperbolic if A is diagonalizable with real eigenvalues; 
we will henceforth assume this to be the case. Let and for 1 < p < m denote the eigenvalues and the 
corresponding right eigenvectors of A with the eigenvalues ordered so that < < . . . < s™. 
Consider the Riemann problem consisting of (2.1) together with initial data 

The solution for i > is piecewise constant with m discontinuities, the pth one proportional to and 
moving at speed . They can be obtained by decomposing the difference ~ qi in terms of the eigenvectors 

qr-qi=Y.aPrP = Y.W^. (2.3) 
P V 

We refer to the vectors W' as waves. Each wave is a jump discontinuity along the ray x = s^t. An example 
solution is pictured in Figure 2.1 for m = 3. For brevity, we will sometimes refer to the Riemann problem 
with initial left state qi and initial right state qr as the Riemann problem with initial states {qi,qr). 

In a finite volume method, it is useful to define notation for the net effect of all left- or right-going waves: 



yt^Ag = ^(sf)- (2.4a) 
p=i 

m 

A+Aq = J2 (s^)"^ (2-4b) 
p=i 
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Here and throughout, {x)^ denotes the positive or negative part of x: 



(x) = min(x, 0) (a;)+ = niax(x, 0). 

The symbols A^Aq, referred to as fluctuations, should be interpreted as single entities that represent the 
net effect of all waves travelling to the right or left. The notation is motivated by the constant coefficient 
linear system (2.1), in which case A^Aq = A^{qr — <?/), where A~ (respectively A'^) is the matrix obtained 
by setting all postive (respectively negative) eigenvalues of A to zero. See [19] or [21] for more details. 

The notation for waves and fluctuations defined in (2.3) and (2.4) can also be used to describe numerical 
solutions of Riemann problems for nonlinear systems if the numerical solver approximates the solution by a 
series of propagating jump discontinuities, which is very often the case. Because the approximate Riemann 
solution for a nonlinear system depends not only on the difference Qr — qi but on the values of the states, 
we will sometimes employ for clarity the notation yVP{qi,qr) to denote the pth wave in the solution of the 
Riemann problem with initial states {qi,qr). In this case the vectors used for the decomposition (2.3) 
are typically eigenvectors of an averaged Jacobian matrix A, and the are the corresponding wave speeds. 
Then A^Aq = A^{qr — qi)- But there are other possible ways to define the r^, and hence the fluctuations. 
For example, using the HLL approximate Riemann solver ]10] would always use only two waves regardless 
of the size of the system. Or, for the case of spatially varying flux functions, the left-going waves may be 
defined by eigenvectors of the Jacobian f'{qi) while the right-going waves are defined by eigenvectors of 
the Jacobian f'{qr)- One could combine these to create a matrix A having this set of eigenvectors and the 
corresponding eigenvalues, but this is not required for implementation of the method. For these reasons we 
use the more general notation A^Aq for the fluctuations deflned by (2.4). 

2.2. First-order Godunov's method. Consider the constant-coefficient linear system in one dimen- 
sion (2.1). Taking a flnite volume approach, we deflne the cell averages (i.e. the solution variables) 

1 r-^h 

Qt{t) = ^ / q{x,t)dx, 
Ax J^ _^ 

* 2 

where the index i and the quantity Ax denote the cell index and the cell size, respectively. To solve the 
linear system (2.1), we initially approximate the solution q{x,t) by these cell averages; that is, at t = to we 
deflne the piecewise-constant function 

q{x,to) = for x e {x^_i,x^^i). (2.5) 

The linear system (2.1) with initial data q consists locally of a series of Riemann problems, one at each 
interface x.^_i. The Riemann problem at Xj„i consists of (2.1) with the piecewise constant initial data 



q{x,0) 



Qj-i X < x^^i 



As discussed above, the solution of the Riemann problem is expressed as a set of waves obtained by decom- 
posing the jump in Q in terms of the eigenvectors of A: 



p 



Let q{x, to + At) denote the exact evolution of q after a time increment At. If we take At small enough 
that the waves from adjacent interfaces do not pass through more than one cell, then we can integrate (2.1) 
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Fig. 2.2. Time evolution oj the reconstructed solution q in cell i. 



over [x^_i^ ^i+i] X [0; At] and divide by Ax to obtain 



1 



(2.7) 



Here should be understood in the sense of distributions. 

We can split the integral above into three parts, representing the Riemann fans from the two interfaces, 
and the remaining piece: 

Aq^dx = / Aq^dx 

X . \ J X . I 



x.,i+s'^At 



Aq^dx + / Aq^dx. (2. 

X. i+s'^At 



The relevant regions are depicted in Figure 2.2. Here we have defined = min(s-'^ i,0) and = 
max(s™ 1 , 0). The third integral in (2.8) vanishes because q(x, At) is constant outside the Riemann fans, by 

* 2 

the definition (2.5). Hence (2.8) reduces to 



'^^Aq^dx^Atj:{sl^) wi^+Atj:{s^^^-wi^ 

-1 p=l - - - 

= At (^+AQ,„i +^-Ag,+ . 



(2.9a) 
(2.9b) 



where the fluctuations A AQ^_^i and A'^AQ^_i are defined by 

m _ 

A-AQ^^.^Y.{<+i) Wf+i' Wf+i (2.10a) 
^+AQ,_i EE^(sP_,) Wf_,, (2.10b) 

Note again that the fluctuations A^AQ^_i and A^AQ^^l are motivated by the idea of a matrix- vector 
product but should be interpreted as single entities that represent the net effect of all waves travelling to the 
right or left. The upper-case Q in the fluctuations is meant to emphasize that they are based on differences 
of cell averages. For instance, the fluctuation A'^AQ^_i corresponds to the effect of right-going waves from 
the Riemann problem with initial states {Qi~i,Qi). 

Substituting (2.9b) into (2.7), we obtain the scheme 

Dividing by At and taking the limit as At approaches zero, we obtain the semi-discrete wave-propagation 
form of the (flrst-order) Godunov's scheme 

^Q^ 1 



dt Ax 



A+AQ,_.+A-AQ,+ .). (2.11) 



Equation (2.11) constitutes a linear system of ordinary differential equations (ODEs) that may be integrated, 
for instance, with a Runge-Kutta method. It is clear from the derivation that this scheme reduces to the 
corresponding flux-differencing scheme when applied to systems written in conservation form. e.g. system 
(1.2). 

2.3. Extension to higher order. The method of the previous section is only flrst-order accurate in 
space. In order to improve the spatial accuracy, we replace the piecewise-constant approximation (2.5) by a 
piecewise-polynomial approximation that is accurate to order p in regions where the solution is smooth: 

q{x,to) = qi{x) for x e {x^_i,x^^i), (2.12) 

where 

q,{x)=q{x,to) + 0{AxP+^). 

Integration of ^4^3; over [xj^_i,Xj^^i] again yields (2.8), but now the third integral is non-zero in general, 
since q is not constant outside the Riemann fans. Deflne 

q^_^ = lim q,ix) i = lim q,{x), (2.13) 

*■ 2 X—^X ^ X >-X~ 

where superscripts L and R refer respectively to the left and the right state of the interface considered. 
Then in place of (2.9), we now obtain 

/ Aq.,dx = At V (si , ) Wf_ , + Ai V (s^ A Wf -t- / Aq^dx (2.14a) 

^At[A+Aq,_.+A-Aq,+ .)+A{q^^,-ql,). (2.14b) 
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The resulting fully-discrete scheme is thus 



We use the notation A^Aq instead of A^AQ because the states in the Riemann problems are not the cell 
averages, but rather the reconstructed interface values. In other words, the fluctuations at x^_i are defined 

by 

rn _|_ 

For instance, the fluctuation A'^Aq^^i corresponds to the effect of right-going waves from the Riemann 
problem with initial states (g^ i , q^_i)- Moreover, we can view the term A{q^ ^ — q^_i) as the sum of both 
the left- and right-going fluctuations resulting from a Riemann problem with initial states (9^1,9^1)- It 
is natural to denote this term, which we refer to as a total fluctuation, by AAqi: 

AAq, = Y: W^l'zf-i^'zti)- 
Dividing by At and taking the limit as At approaches zero, we obtain the semi-discrete scheme 



dt Ax 



A+Aq^_i+A-Aq^^i+AAq,]. (2.15) 



2.4. Spatially varying linear systems. Next we generalize the method to solve one-dimensional 
spatially varying nonconservative linear systems: 

qt+A{x)q^^Q. (2.16) 

We again assume that A is a constant function of x within each cell, so we can write A{x) = Ai{q). In the 

special case that A is the Jacobian matrix of some function /, (2.16) corresponds to the conservation law 

(1.1). Our method can be applied to the general system (2.16) as long as the physically meaningful solution 

to the Riemann problem can be approximated. This may be nontrivial for a nonconservative nonlinear 

problem, as discussed in many papers such as [3, 4, 5, 6, 16]. We do not go into this here, but assume 

that our numerical approach is to be applied to a problem for which the user has a means of solving the 

Riemann problem in terms of discontinuities (waves) 1 propagating at speeds 1 , and hence can define 

* 2 ^2 
fluctuations. This is the case for any linearized Riemann solver, and often for other approximate solvers. 

Then the scheme is given by 



dt Ax 



A+Aq,_i+A-Aqi^i+ I A^q^dx]. (2.17) 



In general, the integral in (2.17) must be evaluated by quadrature; however, for the conservative system 
(1.1), the integral can be evaluated exactly, and is given by 

Aq^dx = f{q^_^r) - f{ql^). (2.18) 

1 
2 
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If the fluctuations are computed using a Roe solver or some other conservative wave-propagation Riemann 
solver, then the flux difference appearing in (2.18) is equal to the sum of fluctuations from a fictitious 
"internal" Riemann problem for the current cell i, just as in the linear case above: 



fiqt+i) - /(Ci) = + A-Aq, = AAq,. (2.19) 



Specifically, the fiuctuations A Aqi are those resulting from the Riemann problem with initial states 



{q^i,q^,i)- Then we can write (2.17) also as 



dQ^ ^ 1 
dt Ax 



A-Aq^^l+A+Aq.^_i+AAq,). (2.20) 



Note that, for the conservative system (1.1), if a Roe solver or an /-wave solver (see Section 2.6) is used, 
then the fiuctuations are equal to the fiux differences 

A-Aq,_.=f.^_^-f{ql^J (2.21) 
A+Aq,_.^fiql,J-f,_., (2.22) 

where /j_ i is the numerical flux at Xj_ i . Thus (2.20) is equivalent to the traditional flux-differencing method 

f^-^ (/...-/._.). (2.23) 

In particular, the scheme is conservative in this case. 

2.5. Capacity-form differencing. In many applications the system of conservation laws takes the 
form 



'i{x)qt + f{q). = 0, (2.24) 

in one space dimension, or 

<x,y)qt + f{q).+9{q)y^0, (2.25) 

in two dimensions, where k is a given function of space and is usually indicated as capacity function (see 
[19]). Systems like (2.24) and (2.25) arise naturally in the derivation of a conservation law, where the flux 
of a quantity is naturally deflned in terms of one variable q, whereas it is a different quantity Kq that is 
conserved. For instance, for the flow in a porous media, k would be the porosity. Note that a capacity 
function can also appear in systems that are not in conservation form, e.g. the quasilinear system (1.3). 

Several approaches can be used to reduce such a system to a more familiar conservation law. One natural 
approach is the capacity- form differencing [19], 

'-W - (-4+A,._. +^-A,.,. +AAq.) , (2.26) 

where Ki is the capacity of the ith cell. This is a simple extension of (2.15) or (2.20) which ensures that 
^ KiQi is conserved (except possibly at the boundaries) and yet allows the Riemann solution to be computed 
based on g as in the case k = 1. 
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2.6. /-wave Riemann solvers and well-balancing. For application to conservation laws, it is desir- 
able to ensure that the wave-propagation discretization is conservative. This can easily be accomplished by 
using an /-wave Riemann solver [1]. Use of /-wave solvers is also useful for problems with spatially varying 
flux function, as well as problems involving balance laws near steady state. 

The idea of the /-wave splitting for (1.1) is to decompose the flux difference f{qr) ~ fill) into waves 
rather than the g-difference used in (2.6), i.e. we decompose the flux difference as a linear combination of 
the right eigenvectors of some Jacobian: 

fiqr) - f{qi) ^Y.P'''-'' = E^''(9i,9r). (2.27) 

p p 

The fluctuations are then deflned as 

p:sP<0 p:sP>0 

Note that the total fluctuation in cell i is given simply by 

-4A,. = /(,^,.)-/(e.). 

An advantage of particular interest is the possibility to include source terms directly into the /-wave 
decomposition. In fact, for balance laws that include non-hyperbolic terms, 

qt + f{q)x = i^iq,^), 

one can easily extend this algorithm by first discretizing the source term to obtain values at the cell 

interfaces and then compute the waves i by splitting 

' 2 

f{qr) - f{qi) - Ax^{qi,qr,x) = ^ /J^'rf = ^ (g,, q,). (2.28) 

p p 

Here ^{qi, qr, x) is some suitable average of ilj{q, x), between the neighboring states. In Bale et al. [1], it has 
been shown that for the second-order FV wave-propagation scheme implemented in Clawpack, the /-wave 
approach is very useful for handling source terms, especially in cases where the solution is close to a steady 
state because it leads to a well-balanced scheme. However, for our high order wave-propagation scheme, 
application of the /-wave algorithm with component-wise or characteristic-wise reconstruction [25] (which 
take no account of the source term) does not lead to a method that is well-balanced, even though the source 
term is accounted for in the Riemann solves which compute A^Aq^_i. This is because with the aforementined 

WENO approaches the reconstructed solution within each cell is not constant (i.e. q^_^i ¥^ qfLi) ^^d 
AAqi ^ 0. In order to obtain a high-order well-balanced scheme, the /-wave Riemann solver is used in 
combination with a /-wave-slope reconstruction (see Section 2.7.3) and a proper Gauss quadrature approach 
for the calculation of the source term. The reconstruction methods considered in this work are presented in 
the next section. 

2.7. Reconstruction. The reconstruction (2.12) should be performed in a manner that yields high 
order accuracy but avoids spurious oscillations near discontinuities. For this purpose, we use weighted 
essentially non-oscillatory (WENO) reconstruction [27]. The spatial accuracy of the method will in general 
be equal to that of the reconstruction. In the present work we employ fifth-order WENO reconstruction. 
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2.7.1. Scalar WENO Reconstruction. WENO reconstruction formulas are typically written in 
terms of the divided differences AQj_|_i/Ax. It is possible to rewrite them in terms of the difference ra- 
tios 

e,_i , ^ ' ' ■ (2.29) 

as long as A(5j_i ^ 0. Then the reconstructed interface values in cell i are given by 

<zf_i =Q^-m-i,2-fc,---,^.-i,fc-i)Ag,_i (2.30a) 
=Q» + '^(^«+i,i-fe,---,^.+ i,fc-2)AQ,_i, (2.30b) 

where is a particular nonlinear function that we will not write out here. The usefulness of (2.30) is that 
it allows WENO reconstruction to be applied to waves directly by redefining 0, as we do below. In the case 
that AQ,_i « (to near machine precision), the value of (f) is set to zero. 

For systems of equations, the simplest approach to reconstruction is component-wise reconstruction, 
which consists of applying the scalar reconstruction approach (2.29)-(2.30) to each element of q. A more so- 
phisticated approach is characteristic-wise reconstruction, in which an eigendecomposition of q is performed, 
followed by reconstruction of each eigencomponent. For problems with spatially-varying coefficients, even 
the characteristic-wise reconstruction may not be satisfying, since it involves comparing coefficients of eigen- 
vectors whose direction in state space varies from one cell to the next. In Clawpack, an alternative kind of 
TVD limiting known as wave limiting has been implemented and shown to be effective for such problems. 

2.7.2. Wave-slope reconstruction. In order to implement a wave-type WENO limiter, we first solve 
a Riemann problem at each interface using the adjacent cell average values Qi-i,Qi as left and right 
states. This results in a set of waves W^ i , which are used solely for the purpose of reconstruction. This 

* 2 

reconstruction is performed by replacing (2.29) by 

^Li. = J7 ' ...n (2.31) 



and replacing (2.30) by 



+ E ^.-1 ■ ■ • ' (^^.,,-2)^U (2.32a) 

p 

e^^Q^ -E^-^(^ri,2-.'---'^ri,.-i)wri. (2-32b) 

This approach takes into account the smoothness of the pth characteristic component of the solution by 
using the information arising from the /c-cell stencils. It is intended to be similar to that used in Clawpack 
[21] and can be conveniently implemented using the same Riemann solvers supplied with Clawpack. 

2.7.3. /-wave-slope reconstruction. Wave-slope reconstruction can also be performed using an /- 
wave Riemann solver. This is useful for computing near-equilibrium solutions of balance laws. The procedure 
is identical to that above, except that the Riemann problem is solved with the /-wave Riemann solver at each 
interface x^_i, using the adjacent cell average values Qi-i,Qi as left and right states. Since the resulting 
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/-waves have the form of a g increment multipHed by the wave speed [1], the waves Z are normahzed by the 
corresponding wave speed before being used for reconstruction: 



WP_i =ZP_i/5^'_i. (2.33) 

* 2 * 2 * 2 

Particular attention must be given to the special situations of = which requires — when the wave- 
propagation approach leads to a conservative algorithm [1]. The reconstruction procedure (2.31)-(2.32) is 
then applied to the waves computed by (2.33). 

For systems with source terms that are at a steady state, assuming that (2.28) results in = for all 
waves, then this /-wave-slope reconstruction will yield to a constant reconstructed function in regions where 
source terms and hyperbolic terms are balanced. Then all fluctuations computed in the update step will 
vanish, so the steady state will be preserved exactly. 

2.8. Time integration. The semi-discrete scheme can be integrated in time using any initial-value 
ODE solver. Herein we use the ten-stage fourth-order strong-stability-preserving Runge-Kutta scheme of 
[11]. This method has a large stability region and a large SSP coefficient, allowing use of a relatively large 
CFL number in practical computations. In all numerical examples of the next section, a CFL number of 
2.45 is used. 

To summarize, the full semi-discrete algorithm used in each Runge-Kutta stage is as follows. 

0. (only if using wave-slope reconstruction) Solve the Riemann problem at each interface x^_i using 
the adjacent cell average values Qi-i,Qi as left and right states. 

1. Compute the reconstructed piecewise function q, and in particular the states q^^i , Q^_^i in each cell, 
using component- wise, characteristic-wise, or wave-slope reconstruction. 

2. At each interface compute the fluctuations A'^Aq^_i and A~Aq^_i by solving the Riemann 
problem with initial states {q^_i,q^_i)- 

^2 * 2 

3. Over each cell, compute the integral J Aq^dx. For conservative systems this is just the total fluctu- 
ation AAqi. 

4. Compute dQ/dt using the semi-discrete scheme (2.17). 

Note again that, for conservative systems, the quadrature in step 3 requires nothing more than evaluating 
and differencing the fluxes. 

2.9. Extension to Two Dimensions. In this section, we extend the numerical wave propagation 
method to two dimensions using a simple dimension-by-dimension approach. The method is applicable to 
systems of the form 

qt + A{x, y)q^ + B{x, y)qy = (2.34) 

on uniform Cartesian grids. 

The 2D analog of the semi-discrete scheme (2.20) is 

- ^^^^^ 

For the method to be high order accurate, the fluctuation terms like A^ lS.q^j_i j should involve integrals over 
cell edges, while the total fluctuation terms like AAqij should involve integrals over cell areas. This can be 
achieved by forming a genuinely multidimensional reconstruction of q and using, e.g.. Gauss quadrature. An 
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implementation following this approach exists in the SharpClaw software. For nonlinear problems containing 
shocks, the genuinely multidimensional reconstruction has been found to be inefficient (at least for some 
simple test problems), as it typically yields only a small improvement in accuracy over the dimension- 
by-dimension scheme given below (which formally only second-order accuarate), but has a much greater 
computational cost on the same mesh. Similar observations have been reported in Zhang et al. [29], where 
both approaches have been throughly tested and compared for linear and nonlinear systems. For problems 
with shocks, at least for the simple test problems presented, the two schemes give comparable resolution on 
the same meshes, despite their difference in formal order of accuracy. 

We now describe the dimension-by-dimension scheme for a single Runge-Kutta stage. We first reconstruct 
piecewise-polynomial functions qj{x) along each row of the grid and qi{y) along each column, by applying a 
ID reconstruction procedure to each slice. We thus obtain reconstructed values 

qf{x,_i)Kq{x,_i,yj) (2.36a) 

qfix.+ i) ~ g(a;,+ i,2/j) (2.36b) 

gf(y,_i) w q{x„yi_i) (2.36c) 

qUy^+^) ~ q{x^,yi+l) (2.36d) 

for each cell The fluctuation terms in (2.35) are determined by solving Riemann problems between the 
appropriate reconstructed values; for instance B~ Aq^^ jj^i is determined by solving a Riemann problem in 
the y-direction with initial states {q^j^i,q^'j^i)- In the case of conservative systems or piecewise-constant 
coefficients, the total fluctuation terms AAqi_j and BAqij can be similarly determined by summing the left- 
and right-going fluctuations of an appropriate Riemann problem. Thus, for instance, BAqij is determined 
by solving Riemann problem in the j/-direction with initial states {q^ _i,q^-,i)- 

3. Numerical applications. In this section we present results of numerical tests using the wave 
propagation methods just described. The examples included are chosen to emphasize the advantages of the 
wave propagation approach. We make some comparisons with the well-known second-order wave propagation 
code Clawpack [22, 21[. 

3.1. Acoustics. In this section, the high-order wave propagation algorithm is applied to the ID equa- 
tions of linear acoustics in piecewise homogeneous materials: 

Pt + K{x,y){ux+Vy) ^0 (3.1a) 
.. + ^p.=0 (3.1b) 

+ = 0. (3.1c) 

Here p is the pressure and u, v are the x- and y- velocities, respectively. The coefficients p and K, which vary 
in space, are the density and bulk modulus of the medium. We will also refer to the sound speed c — yj Kj p. 
Notice that in general since p varies in space, the last two equations above are not in conservation form. This 
test case demonstrates that the proposed approach is able to solve hyperbolic system of equations written 
in nonconservative form. Of course, this system can be written in conservative form as follows: 

e* - {ux + Wy) = (3.2a) 

p{x,y)ut~{K{x,y)t\ = ^ (3.2b) 

p{x, y)vt - iK{x, y)e)y = 0, (3.2c) 
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Where e — ~p/K is the strain. As we will see, the latter form may be advantageous in terms of the accuracy 
that can be obtained. 

The grid is chosen so that the material is homogeneous in each computational cell, and an exact Riemann 
solver is used at each interface; for details of this solver see e.g. [9]. 

3.1.1. One-dimensional acoustics. We first consider one-dimensional acoustic waves in a piecewise- 
constant medium with a single interface. Namely, we solve (3.1) on the interval x E [—10, 10] with 



{pi,ci) x<0 

{pr, Cr) X > 



We measure the convergence rate of the solution in order to verify the order of accuracy for smooth solutions. 
The initial condition is a compact, six-times differentiable purely right-moving pulse: 

pix, 0) = ii---o)^anix-xo)+ar ^^^ _ 



u{x,0) ^p{x,0)/Z{x) 



where 



^{X - Xq) 



for — xol > a) 

1 for |a; — xol < a. 



with Xf) = —4 and a = 1, and Z{x) ~ K{x)p{x). This condition was chosen to be sufficiently smooth to 
demonstrate the design order of the scheme and to give a solution that is identically zero at the material 
interface at the initial and final times. 

Table 3.1 shows Li errors and convergence rates for propagation in a homogeneous medium with pi = 
ci — Pr — Cr — 1. Here we use componentwise reconstruction. Specifically, we compute 



EL,=AxJ2\Q^~Q^ 



(3.3) 



where Qi is a highly accurate solution cell average computed by characteristics or by using a very fine 
grid. For the acoustics problems in this section, Q is computed using characteristics and adaptive Gauss 
quadrature. Table 3.1 indicates that in each case, the order of convergence is approximately equal to the 
design order of the discretization. 

Table 3.1 
Errors for homogeneous acoustics test 





SharpClaw 


Clawpack 


mx 


Error Order 


Error Order 


200 


3.60e-02 


4.10e-02 


400 


3.65e-03 3.30 


1.30e-02 1.66 


800 


1.85e-04 4.31 


3.61e-03 1.85 


1600 


7.35e-06 4.65 


8.94e-04 2.01 



To test the accuracy in the presence of discontinuous coefficients we take 



P; = C; = 1 Pr=4: Cr = 1/2, 
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with an impedance ratio of Z^jZi ~ 2. As was noted in [1], this system can also be solved in the conservative 
form (3.2) using the /-wave approach. We include results of this approach, where we have also performed 
characteristic-wise rather than component-wise reconstruction. Results are shown in Table 3.2. In this case 
all schemes exhibit a convergence rate below the formal order, even though the initial and final solutions are 
smooth. To investigate this further, we repeat the same test with a wider pulse by taking a = 4. Results are 
shown in Table 3.3. 

For the latter test, we observe a convergence rate of approximately two for SharpClaw, one for Claw- 
pack, and five for SharpClaw using the /-wave approach and characteristic-wise reconstruction. The last 
convergence rate is remarkable, considering that the solution is not differentiable when it passes through 
the material interface. Further investigation of the accuracy of this approach for more complicated prob- 
lems with discontinuous coefficients is ongoing. In tests not shown here, Clawpack achieves approximately 
second-order accuracy when used with an /-wave Riemann solver for this problem. 



Table 3.2 

Errors for acoustics interface with narrow pulse 





SharpClaw 


SharpClaw 


/-wave 


Clawpack 


mx 


Error 


Order 


Error 


Order 


Error Order 


200 


2.10e-01 


9.50e-02 


1.98e-01 


400 


5.98e-02 


1.81 


1.42e-02 


2.74 


7.26e-02 1.45 


800 


1.25e-02 


2.26 


1.42 e-03 


3.32 


2.21e-02 1.71 


1600 


1.17e-03 


3.42 


1.20e-04 


3.56 


7.86e-03 1.49 


L 


^ Errors for 


Table 3.3 
acoustics interface probl 


em with wide pulse (a=4) 




SharpClaw 


SharpClaw 


/-wave 


Clawpack 


mx 


Error 


Order 


Error 


Order 


Error Order 


200 


9.67e-03 


5.01e-03 


5.23e-02 


400 


2.01e-03 


2.27 


4.63e-04 


3.44 


2.32e-02 1.17 


800 


4.89e-04 


2.04 


2.51e-05 


4.36 


1.09e-02 1.09 


1600 


1.22e-04 


2.00 


6.49e-07 


5.12 


5.26e-03 1.05 



3.1.2. A Two-dimensional sonic crystal. In this section we model sound propagation in a sonic 
crystal. A sonic crystal is a periodic structure composed of materials with different sounds speeds and 
impedances. The periodic inhomogeneity can give rise to bandgaps - frequency bands that are completely 
reffected by the crystal. This phenomenon is widely utilized in photonics, but its significance for acoustics has 
only recently been considered. Photonic crystals can be analyzed quite accurately using analytic techniques, 
since they are essentially infinite size structures relative to the wavelength of the waves of interest. In 
contrast, sonic crystals are typically only a few wavelengths in size, so that the effects of their finite size 
cannot be neglected. For more information on sonic crystals, see for instance the review paper [24]. 

We consider a square array of square rods in air with a plane wave disturbance incident parallel to one 
of the axes of symmetry. The array is infinitely wide but only eight periods deep. The lattice spacing is 10 
cm and the rods have a cross-sectional side length of 4 cm, so that the filling fraction is 0.16. This crystal is 
similar to one studied in [26], and it is expected that sound waves in the 1200-1800 Hz range will experience 
severe attenuation in passing through it, while longer wavelengths will not be significantly attenuated. 
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Fig. 3.1. Pressure in the sonic crystal for a long wavelength plane wave incident from the left. 




Fig. 3.2. Pressure in the sonic crystal for a long wavelength plane wave incident from the left. 



A numerical instability very similar to that observed in ID simulations in [8, 9] was observed when the 
standard Clawpack method was applied to this problem. The fifth-order WENO method with characteristic- 
wise limiting showed no such instability. 

Figure 3.1 shows a snapshot of the RMS pressure distribution in space for a plane wave with k — 15 
incident from the left. The RMS (root mean square) pressure is computed as follows: 



PRMs{x,y) = d - J p'^{x,y,t)dt. (3.4) 

This wave has a frequency of about 800 Hz, well below the partial band gap. As expected, the wave passes 
through the crystal without significant attenuation. In Figure 3.2, the pressure is plotted along a slice in the 
cc-direction approximately midway between rows of rods. 

Figure 3.3 shows a snapshot of the RMS pressure distribution in space for an incident plane wave with 
with frequency 1600 Hz, inside the partial bandgap. Notice that the wave is almost entirely refiected, 
resulting in a standing wave in front of the crystal. Figure 3.4 shows the RMS pressure along a slice in the 
x-direction. 

3.2. Nonlinear Elasticity in a Spatially Varying Medium. In this section we consider a more dif- 
ficult test involving nonlinear wave propagation and many material interfaces. This problem was considered 
previously in [17] and studied extensively in [18]. Solitary waves were observed to arise from the interaction 
of nonlinearity and an effective dispersion due to material interfaces in layered media. 

Elastic compression waves in one dimension are governed by the equations 

et{x,t) — Ux{x,t) — (3.5a) 
{p{x)u{x, t))t - cr(e(a;, t),x)x = 0. (3.5b) 
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Fig. 3.3. Snapshot of RMS pressure distribution in space in the sonic crystal for a plane wave incident from the left. 




where e is the strain, u the velocity, p the density, and a the stress. This is a conservation law of the form 
(1.1), with 

Note that the density and the stress-strain relationship vary in x. The Jacobian of the flux function is 

/'<')^ -r' 

In the case of the linear stress-strain relation a{x) = K{x)e{x), (3.5) is equivalent to the one-dimensional 
form of the acoustics equations considered in the previous section. 
We consider the piecewise constant medium studied in [17, 18]: 

(n(T) K(t)) - / ^) if < 2^ < (i + 5) for some integer j 
(p[x),K(x)) - <^ (4,4) otherwise, 

with exponential stress-strain relation 

cr(e, x) = exp{K{x)e) - 1. (3.7) 

The initial condition is uniformly zero, and the boundary condition at the left generates a half-cosine pulse. 

We solve this problem using the /-wave solver developed in [17]. Figure 3.5 shows a comparison of 
results using Clawpack and SharpClaw on this problem, with 24 cells per layer. The SharpClaw results are 
significantly more accurate. 
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strain at time 600.0000 



Stress at time 600.0000 



0.0 nun 




(a) Strain. 



(b) Stress. 



Fig. 3.5. Comparison of Clawpack (red circles) and SharpClaw (blue squares) solution of the stegoton problem using 24 
cells per layer. For clarity, only every third solution point is plotted. The black line represents a very highly resolved solution. 

Solutions of (3.5) are time-reversible in the absence of shocks. As discussed in [12, 14], the effective 
dispersion induced by material inhomogeneities suppresses the formation of shocks, for small amplitutde 
initial and boundary conditions, rendering the solution time-reversible for very long times. This provides a 
useful numerical test. We solve the stegoton problem numerically up to time T, then negate the velocity and 
continue solving to time 2T. The solution at any time 2T — tg, with to < T, should be exactly equal to the 
solution at to. We take T — 600 and to = 60. Figure 3.6(a) shows the solution obtained using SharpClaw 
on a grid with 24 cells per layer. The t = 1140 solution (squares) is in excellent agreement with the < = 60 
solution (solid line). In fact, the maximum point-wise difference has magnitude less than 2 x 10^^. Using 
a grid twice as fine, with 48 cells per layer, reduces the point- wise difference to 1 x 10^"^. The Clawpack 
solution, computed on the same grid (24 cells per layer), is shown in Figure 3.6(b). Again, the SharpClaw 
solution is noticeably more accurate. For a more detailed study of this time-reversibility test, we refer to 
[14]. 

3.3. Shallow Water Flow. In order to show the capabilities of the proposed scheme to deal with 
nonlinear problems with source terms, the shallow water equations are also considered. The conservative 
form of the depth-averaged equations of mass and momentum in two space dimensions can be written as 
follows: 



ht + {hu)x + {hv)y = 
{hu)t + ( ]^hu^ + \9^j + {huv)y = -ghh^ 



{hv)t + {huv).^ + ( ^hv^ + ]^gh? 



-ghby 



(3.8a) 
(3.8b) 

(3.8c) 



where h, u and v are the depth of the fluid and the velocity components in x and y directions, respectively. 
The function b{x,y) is the bottom elevation and g is constant for gravitational acceleration. 

Two test cases are presented: a radially symmetric dam-break problem over a flat bottom topography 
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stress at time 1140.0000 




Stress at time 1140.0000 




(a) SharpClaw. (b) Clawpack. 

Fig. 3.6. Comparison of forward solution (black line) and time-reversed solution (symbols). 



(h — 0) and a small perturbation of a steady state over a hump. A Roe solver with entropy fix is used in 
both cases. 

3.3.1. Radial dam-break problem. This problem consists in computing the flow induced by the 
instantaneous collapse of an idealized circular dam. It is widely used to benchmark various numerical 
techniques that tend to simulate interfacial flows and impact problems. 

The dam, which is an infinitesimally thin circular wall with a radius of 0.5, is located at the center of a 
square computational domain of 2.5 x 2.5. The water level is initially h = 2.0 inside the circle and h = 1.0 
outside. The initial solution has zero discharge, i.e. both velocity components are zero. The dam is removed 
at time t — Q. This tests the ability of the method to compute the 2D propagation of nonlinear waves and 
the extent to which symmetry is preserved in the numerical solution. In the presence of radial symmetry, 
system (3.8) can be recast in the following form: 



ht + {hU)r = 



{hU)t + ( \hU^ + \gh^ 



hU 
r 



(3.9a) 
(3.9b) 



where h is still the depth of the fluid, whereas U and r are the radial velocity and the radial position. 
An important feature of these equations is the presence of a source term, which physically arises from the 
fact that the fluid is spreading out and it is impossible to have constant depth and constant non-zero radial 
velocity. 

A first comparison between SharpClaw and Clawpack is performed by solving the ID system (3.9) on 
the interval < r < 2.5. A wall boundary condition and non-reflecting boundary condition are imposed at 
the left and the right boundaries, respectively. The final time for the analysis is taken to be i = 1. The 
classical g-wave Riemann solver based on the Roe linearization is used to solve the Riemann problem at each 
interface (see for instance [21] for details), where the left and the right states are computed by using the 
characteristic-wise WENO reconstruction. The gravitational acceleration is set to g = 1. A highly resolved 
solution obtained with Clawpack on a grids with 25, 600 cells is used as a reference solution. 
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Pointwise error norm at time 1 .0000 
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Fig. 3.7. Pointwise absolute error for the water height on a grid with 800 cells. SharpClaw solution: continuous line; 
Clawpack solution: dashed line. 



It is well-known that formal order of accuracy is lost in a shock wave propagating in a coupled system 
of equation [23] and in general it is reduced to first-order. However, if we plot the difference between the 
computed solution available at the cell's center and the reference solution conservatively averaged on the 
same grid, i.e. Ei = \Qi — Qi\, then we can visualize where the errors are largest as well as their spatial 
structure. Figure 3.7 shows this difference for the water height (h) on a grid with 800 cells. The reference 
solution at i = 1 is shown by the solid line in Figure 3.8. The largest errors in both solutions are near the 
shocks. In the smooth regions, the SharpClaw solution is more accurate than that of the Clawpack code. 

Next we consider the same problem using the full 2D equations (3.8). The SharpClaw and Clawpack 
codes are tested on two grids with 125 x 125 and 500 x 500 cells. The final time for the analysis is again 
taken to be i = 1. 

Figure 3.8 shows the water height h computed with SharpClaw at each cell's center and t = 1.0 as a 
function of the radial position. The radius is measured respect to the center of the initial condition. The 
ID reference solution used before is also plotted for comparison. It can be seen that the scheme preserves 
a good radial symmetry, though it cannot resolve the shock near the origin. The grid is in fact too coarse. 
Clawpack results are not shown in this figure but indicate similar accuracy and similarly good symmetry. 

The solutions obtained on the finer grid (500 x 500 cells) are shown in Figure 3.9. The effect of the grid 
refinement is clearly visible. In fact, the solutions gets close to the reference solution. However, the density 
of the grid near the origin is still coarse to resolve accurately the shock near the origin. 

3.3.2. Perturbation of a steady state solution. Conservation laws with source terms often have 
steady states in which the flux gradient are non-zero but exactly balanced by source terms. A good numerical 
scheme, should be able to preserve such steady states and calculate accurately small perturbations around 
these conditions. A classical benchmark test case to investigate these properties is the small perturbation of 
a 2D steady state water given by LeVeque [20] . 

System (3.8) is solved in a rectangular domain [0, 2] x [0, 1], with a bottom topography characterized by 
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Fig. 3.8. Solution for the 2D radial dam-break problem on a grid with 125 x 125 cells, plotted as a function of radius. 
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Fig. 3.9. Solution for the 2D radial dam-break problem on a grid with 500 x 500 cells, plotted as a function of radius. 

an isolated elliptical shaped hump: 

b{x,y) = 0.8exp(-5(a; - 0.9)^ - 50(y - 0.5)^). 

The surface is initially flat with h{x,y,0) = 1 — b{x,y) except for 0.05 < a; < 0.15, where h is perturbed 
upward by e = 0.01. The initial discharge in both direction is zero, i.e. hu{x,y,0) = hv{x,y,0) = 0. Non- 
reflecting (i.e., zero-extrapolation) conditions are imposed at all boundaries. The gravitational acceleration 
is set to 9.81. 

An effort was made to achieve a well-balanced scheme using the /-wave approach combined with 
component-wise or characteristic-wise WENO reconstruction, but this was unsuccessful. This is not sur- 
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Surface level at time t = 0.06000000 ,„ Surface level at time t = 0.12000000 




Fig. 3.10. Contour of the surface level h + b at time t = 0.06 and t = 0.12. Component-wise reconstruction approach. 
Contour levels: 0.99942 : 0.000238 : 1.00656. 

prising, since the algorithm begins by reconstructing a non-constant function. Figure 3.10 shows the contour 
levels of the solution at i = 0.06 and t = 0.12 on a fine grid with 600 x 300 cells, obtained with the /-wave 
Riemann solver and the component- wise reconstruction approach as a building block for the WENO scheme. 
The scheme is not well-balanced and spurious waves are generated around the hump. Similar results are 
obtained using characteristic-wise reconstruction. 

In order to balance the scheme, the /-wave-slope reconstruction introduced in Section 2.7 is used instead. 
In this approach, the WENO reconstruction is applied to waves computed by solving Riemann problems with 
the /-wave solver. When the source term is included in these Riemann problems, the resulting waves vanish 
exactly. Figure 3.11 shows the surface level a cross section along y — 0.5 at time t = 0.06 computed with 
both reconstruction approaches (and the /-wave Riemann solver) on a uniform mesh with 600 x 300. It 
is seen that the /-wave-slope reconstruction method keeps the surface fiat, whereas the component- wise 
reconstruction introduces spurious waves which have an amplitude of the order of the disturbance that we 
want to resolve. 

Figure 3.12 shows the solution on two uniform meshes with 200 x 100 cells and 600 x 300 cells, computed 
using the /-wave-slope reconstruction approach. The results clearly indicate that the detailed structure of 
the evolution of such a small perturbation is resolved well even with the relatively coarse mesh. In addition, 
these results agree with those reported in [20]. 

4. Conclusions. We have presented a general approach to extending the finite volume wave propaga- 
tion algorithm to high order accuracy in one and two dimensions. The algorithm is based on a method- 
of-lines approach, wherein the semi-discrete scheme relies on high order reconstruction and computation 
of fiuctuations, including a total fluctuation term arising inside each cell. By using WENO reconstruction 
and strong stability preserving time integration, high order accurate non-oscillatory results are obtained, as 
demonstrated through a variety of test problems. 

This algorithm has several desirable features. Like the second-order wave propagation algorithms in 
Clawpack [19], it is applicable to hyperbolic PDEs including linear nonconservative systems and nonlinear 
systems with spatially varying fiux function. It has been shown to achieve high order accuracy even for 
problems with discontinuous coefficients. Finally, the algorithm can be adapted to give a well-balanced 
scheme for balance laws by use of the /-wave approach and a new wave-slope reconstruction technique. 

Hyperbolic systems of equations with both smooth and non-smooth solution have been used to test 
the properties and the capabilities of the proposed method. The schemes have been compared for linear 
acoustics and nonlinear elasticity problems in heterogeneous media and for the shallow water equations 
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Surface level at time t = 0.06000000, cross section along y=0.5 
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Fig. 3.11. Surface level h + b along a cross section at y = 0.5 and time t = 0.06. Solid line: f -wave Riemann solver and 
component-wise reconstruction; dashed line: f-wave Riemann solver and f -wave-slope reconstruction. 

with and without bottom topography. Two types of Riemann solver have been used, i.e the classical (q-) 
wave algorithm and the /-wave approach. The new scheme performed well for all the test cases. It gives 
significantly better accuracy than Clawpack (on the same grid) for smooth problems. 

In two dimensions, the presented dimension-by-dimension reconstruction approach is formally only 
second-order accurate (see [29]); however, it gives improved accuracy over the second-order scheme im- 
plemented in Clawpack for the test problems considered. Further investigation of different approaches to 
multidimensional reconstruction for problems containing both shocks and rich smooth flow structures is a 
topic of future research. 
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